Plots and analyses produced for The Indus Village model (Quaternary 2021)
Andreas Angourakis
29/03/2021
This file and all other referenced in the code can be found at the repository: https://github.com/Andros-Spica/Quaternary_Angourakis-et-al-2021
Load packages:
require(reshape2)## Loading required package: reshape2require(ggplot2)## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.6.3require(patchwork)## Loading required package: patchwork## Warning: package 'patchwork' was built under R version 3.6.3require(colorspace)## Loading required package: colorspace## Warning: package 'colorspace' was built under R version 3.6.3require(grid)## Loading required package: gridrequire(gridExtra)## Loading required package: gridExtra## Warning: package 'gridExtra' was built under R version 3.6.3require(png)## Loading required package: pngDefine lookup functions to access data between tables using identifiers or ‘keys’:
getRelatedVariableWithKey <- function(keyInCallingTable, keyInCalledTable, variableInCalledTable)
{
return(
sapply(keyInCallingTable, function(x) variableInCalledTable[match(x, keyInCalledTable)])
)
}
getRelatedVariableWithTwoKeys <- function(firstKeyInCallingTable, secondKeyInCallingTable,
firstKeyInCalledTable, secondKeyInCalledTable, variableInCalledTable)
{
temp <- 1:length(firstKeyInCallingTable)
for (i in 1:length(temp))
{
temp[i] <- variableInCalledTable[firstKeyInCalledTable == firstKeyInCallingTable[i] &
secondKeyInCalledTable == secondKeyInCallingTable[i]]
}
return(temp)
}Define the format of image files to be generated (only “png” and “eps”):
listOfImageFormats <- c("eps", "png")Crop table
Load crop table:
cropTable <- read.csv("models/cropsTable.csv", skip = 4)Clean empty columns (generated with automatic names starting with X):
cropTable <- cropTable[,names(cropTable)[grep("^(?!X).*$", names(cropTable), perl = TRUE)]]Reorder crops for a more logical display:
WARNING: doing this manually makes the script sensitive to the value of crop-selection in NetLogo. Thus, this chunk must be updated if crop-selection is changed before running experiments.
cropTable$crop <- factor(cropTable$crop, levels = c("proso millet", "pearl millet", "rice", "barley", "wheat 1", "wheat 2"))Get vector of colours to represent crops:
cropColours <- qualitative_hcl(nlevels(cropTable$crop))
# alternative using only base R graphics
#cropColours <- rainbow(nlevels(cropTable$crop), s = 0.8, v = 0.8, end = 0.85)knitr::kable(cropTable)| crop | species | cultivar | T_sum | HI | I_50A | I_50B | T_base | T_opt | RUE | I_50maxH | I_50maxW | T_heat | T_ext | S_CO2 | S_water | sugSowingDay | sugHarvestingDay | root.zone.depth..mm. |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| wheat 1 | Triticum aestivum | Yecora Rojo | 2200 | 0.36 | 480.0 | 200.00 | 0 | 15 | 1.24 | 100 | 25 | 34 | 45 | 0.08 | 0.40 | 330 | 105 | 1000 |
| wheat 2 | Triticum aestivum | Batten | 2150 | 0.34 | 280.0 | 50.00 | 0 | 15 | 1.24 | 100 | 25 | 34 | 45 | 0.08 | 0.40 | 330 | 105 | 1000 |
| rice | Oryza sativa | IR72 | 2300 | 0.47 | 850.0 | 200.00 | 9 | 26 | 1.24 | 100 | 10 | 34 | 50 | 0.08 | 1.00 | 170 | 320 | 400 |
| barley | Hordeum vulgare | hypothetical; based on DSSAT v4.7 | 1762 | 0.42 | 350.0 | 170.00 | 0 | 15 | 1.24 | 100 | 20 | 34 | 45 | 0.08 | 0.40 | 330 | 100 | 1000 |
| pearl millet | Pennisetum glaucum | hypothetical; based on DSSAT v4.7 - MLCER047.CUL - Middle Variety | 1220 | 0.25 | 245.0 | 120.00 | 10 | 33 | 1.90 | 100 | 5 | 35 | 47 | 0.08 | 0.05 | 200 | 320 | 1000 |
| proso millet | Panicum miliaceum L. | hypothetical | 1328 | 0.29 | 157.3 | 96.75 | 0 | 18 | 3.00 | 100 | 130 | 34 | 45 | 0.08 | 0.05 | 107 | 220 | 1000 |
Experiment 1
Experiment 1 simulates the growth and yield of crops included in ‘cropTable’ using the NetLogo implementation of SIMPLE crop model (Zhao et al. 2019) adapted as a module of the Indus Village model under the default parameter configuration, which is aimed to approximate the conditions in Haryana, NW India.
Loading and preparation
Load data:
yieldData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_type-of-experiment=user-defined_experiment-name=exp1_initRandomSeed=0.csv")Convert yield = 0 to NA when meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):
yieldData$yield[is.na(yieldData$meanARID_grow)] <- NAForce the same order of crops used by cropTable:
yieldData$crop <- factor(yieldData$crop, levels = levels(cropTable$crop))Calculate the range of ARID and productivity (yield):
minARID = round(min(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
maxARID = round(max(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
minYield = round(min(yieldData$yield, na.rm = TRUE), digits = -1)
maxYield = round(max(yieldData$yield, na.rm = TRUE), digits = -1)Load daily simulated weather and ARID data:
weatherData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_weather_type-of-experiment=user-defined_experiment-name=exp1_initRandomSeed=0.csv")As the reference real-world data, we use the data downloaded at NASA´s POWER access viewer (power.larc.nasa.gov/data-access-viewer/) selecting the user community ‘Agroclimatology’ between 01/01/1984 and 31/12/2007. The exact location used as a reference is:
- Rakhigarhi, Haryana, India (Latitude: 29.1687, Longitude: 76.0687)
We selected the ICASA Format’s parameters:
- Precipitation (PRECTOT)
- Wind speed at 2m (WS2M)
- Relative Humidity at 2m (RH2M)
- Dew/frost point at 2m (T2MDEW)
- Maximum temperature at 2m (T2M_MAX)
- Minimum temperature at 2m (T2M_MIN)
- All sky insolation incident on a horizontal surface (ALLSKY_SFC_SW_DWN)
- Temperature at 2m (T2M)
Load Soil Water (or Soil Water Balance) model code to calculate ARID given real-world data:
source("FigWeather_input/soilWaterModel/estimateETr.R")
source("FigWeather_input/soilWaterModel/watbal.model.R")Load and format the real-world dataset so that they can be used in the Soil Water model:
weatherData_rakhigarhi <- watbal.weather.file("FigWeather_input/POWER_SinglePoint_Daily_19840101_20071231_029d17N_076d70E_5b401917.csv",
year = 1984:2007)Estimate reference evapotranspiration using FAO recomendations (Penman-Monteith method):
weatherData_rakhigarhi$ETr <- estimateETr(weatherData_rakhigarhi$I,
weatherData_rakhigarhi$T2M,
weatherData_rakhigarhi$Tmax,
weatherData_rakhigarhi$Tmin,
weatherData_rakhigarhi$T2MDEW,
weatherData_rakhigarhi$WS2M)Calculate ARID in real-world data:
Use Soil Water model to calculate soil water content proportion or WATp (mm mm-1 day-1) and ARID coefficient.
The first step is to set all soil parameters, here assumed to be constant (the same values used in experiment 1 simulations):
watbalParameters <- watbal.define.param()Run the model inputing parameters and dataset:
watbalOutput_rakhigarhi <- watbal.model(watbalParameters,
weatherData_rakhigarhi,
typeOfParValue = "nominal")watbal.modeliterates for every day in the dataset calculating soil water content, absolute and proportion (WAT, WATp), and the respective ARID coefficient.
Show parameter values used here in Soil Water model (both in watbal.model and in the SIMPLE crop model). In experiemnt 1, they should be the same:
knitr::kable(
cbind(
parameters = c("elevation", "albedo", "CN", "DC", "z", "FC", "WHC", "WP", "MUF"),
SIMPLE = as.vector(unlist(yieldData[1, c("elevation", "albedo", "CN", "DC", "z", "FC", "WHC", "WP", "MUF")])),
watbal = c(200, # these values of elevation and albedo are the default in estimateETr
0.23,
as.vector(watbalParameters[1, c("CN", "DC", "z", "FC", "WHC", "WP", "MUF")]))
)
)| parameters | SIMPLE | watbal |
|---|---|---|
| elevation | 200 | 200 |
| albedo | 0.23 | 0.23 |
| CN | 65 | 65 |
| DC | 0.55 | 0.55 |
| z | 400 | 400 |
| FC | 0.21 | 0.21 |
| WHC | 0.15 | 0.15 |
| WP | 0.06 | 0.06 |
| MUF | 0.096 | 0.096 |
Summaries per day of year
Create data frames containing the daily summary statistics of simulations and real-world data:
- Create names vector (used for both datasets)
weatherSummaryNames <- c("dayOfYear",
"solarRadiation.mean",
"solarRadiation.sd",
"solarRadiation.max",
"solarRadiation.min",
"solarRadiation.error",
"temperature.mean",
"temperature.sd",
"temperature.max",
"temperature.min",
"temperature.error",
"maxTemperature.mean",
"maxTemperature.max",
"maxTemperature.min",
"maxTemperature.error",
"minTemperature.mean",
"minTemperature.max",
"minTemperature.min",
"minTemperature.error",
"temperature.lowerDeviation",
"temperature.lowerDeviation.error",
"temperature.upperDeviation",
"temperature.upperDeviation.error",
"precipitation.mean",
"precipitation.max",
"precipitation.min",
"precipitation.error",
"ARID.mean",
"ARID.max",
"ARID.min",
"ARID.error")- Simulation data
weatherSummary <- vector("list", length(weatherSummaryNames))
names(weatherSummary) <- weatherSummaryNames
# OBS: the lines above produce an ERROR related to names that is inconsequential
for (day in 1:365)
{
weatherSummary$dayOfYear <- c(weatherSummary$dayOfYear, day)
tempData_thisDay <- weatherData[weatherData$currentDayOfYear == day,]
# solar radiation
weatherSummary$solarRadiation.mean <- c(
weatherSummary$solarRadiation.mean,
mean(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.sd <- c(
weatherSummary$solarRadiation.sd,
sd(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.max <- c(
weatherSummary$solarRadiation.max,
max(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.min <- c(
weatherSummary$solarRadiation.min,
min(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.error <- c(
weatherSummary$solarRadiation.error,
qt(0.975,
length(tempData_thisDay$solarRadiation) - 1) *
sd(tempData_thisDay$solarRadiation, na.rm = TRUE) /
sqrt(length(tempData_thisDay$solarRadiation)))
# temperature
## daily mean
weatherSummary$temperature.mean <- c(
weatherSummary$temperature.mean,
mean(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.sd <- c(
weatherSummary$temperature.sd,
sd(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.max <- c(
weatherSummary$temperature.max,
max(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.min <- c(
weatherSummary$temperature.min,
min(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.error <- c(
weatherSummary$temperature.error,
qt(0.975,
length(tempData_thisDay$temperature) - 1) *
sd(tempData_thisDay$temperature, na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature)))
## daily max
weatherSummary$maxTemperature.mean <- c(
weatherSummary$maxTemperature.mean,
mean(tempData_thisDay$temperature_max, na.rm = TRUE))
weatherSummary$maxTemperature.max <- c(
weatherSummary$maxTemperature.max,
max(tempData_thisDay$temperature_max, na.rm = TRUE))
weatherSummary$maxTemperature.min <- c(
weatherSummary$maxTemperature.min,
min(tempData_thisDay$temperature_max, na.rm = TRUE))
weatherSummary$maxTemperature.error <- c(
weatherSummary$maxTemperature.error,
qt(0.975,
length(tempData_thisDay$temperature_max) - 1) *
sd(tempData_thisDay$temperature_max, na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_max)))
## daily min
weatherSummary$minTemperature.mean <- c(
weatherSummary$minTemperature.mean,
mean(tempData_thisDay$temperature_min, na.rm = TRUE))
weatherSummary$minTemperature.max <- c(
weatherSummary$minTemperature.max,
max(tempData_thisDay$temperature_min, na.rm = TRUE))
weatherSummary$minTemperature.min <- c(
weatherSummary$minTemperature.min,
min(tempData_thisDay$temperature_min, na.rm = TRUE))
weatherSummary$minTemperature.error <- c(
weatherSummary$minTemperature.error,
qt(0.975,
length(tempData_thisDay$temperature_min) - 1) *
sd(tempData_thisDay$temperature_min, na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_min)))
## daily upper and lower deviation
weatherSummary$temperature.lowerDeviation <- c(
weatherSummary$temperature.lowerDeviation,
mean(tempData_thisDay$temperature -
tempData_thisDay$temperature_min))
weatherSummary$temperature.lowerDeviation.error <- c(
weatherSummary$temperature.lowerDeviation.error,
qt(0.975,
length(tempData_thisDay$temperature_min) - 1) *
sd(tempData_thisDay$temperature -
tempData_thisDay$temperature_min,
na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_min)))
weatherSummary$temperature.upperDeviation <- c(
weatherSummary$temperature.upperDeviation,
mean(tempData_thisDay$temperature_max -
tempData_thisDay$temperature))
weatherSummary$temperature.upperDeviation.error <- c(
weatherSummary$temperature.upperDeviation.error,
qt(0.975,
length(tempData_thisDay$temperature_max) - 1) *
sd(tempData_thisDay$temperature_max -
tempData_thisDay$temperature,
na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_max)))
# precipitation
weatherSummary$precipitation.mean <- c(
weatherSummary$precipitation.mean,
mean(tempData_thisDay$RAIN, na.rm = TRUE))
weatherSummary$precipitation.max <- c(
weatherSummary$precipitation.max,
max(tempData_thisDay$RAIN, na.rm = TRUE))
weatherSummary$precipitation.min <- c(
weatherSummary$precipitation.min,
min(tempData_thisDay$RAIN, na.rm = TRUE))
weatherSummary$precipitation.error <- c(
weatherSummary$precipitation.error,
qt(0.975,
length(tempData_thisDay$RAIN) - 1) *
sd(tempData_thisDay$RAIN, na.rm = TRUE) /
sqrt(length(tempData_thisDay$RAIN)))
# ARID
weatherSummary$ARID.mean <- c(
weatherSummary$ARID.mean,
mean(tempData_thisDay$ARID, na.rm = TRUE))
weatherSummary$ARID.max <- c(
weatherSummary$ARID.max,
max(tempData_thisDay$ARID, na.rm = TRUE))
weatherSummary$ARID.min <- c(
weatherSummary$ARID.min,
min(tempData_thisDay$ARID, na.rm = TRUE))
weatherSummary$ARID.error <- c(
weatherSummary$ARID.error,
qt(0.975,
length(tempData_thisDay$ARID) - 1) *
sd(tempData_thisDay$ARID, na.rm = TRUE) /
sqrt(length(tempData_thisDay$ARID)))
}
rm(day, tempData_thisDay)
weatherSummary <- data.frame(weatherSummary)- Real-world data
weatherSummary_rakhigarhi <- vector("list", length(weatherSummaryNames))
names(weatherSummary_rakhigarhi) <- weatherSummaryNames
# OBS: the lines above produce an ERROR related to names that is inconsequential
for (day in 1:366)
{
weatherSummary_rakhigarhi$dayOfYear <- c(weatherSummary_rakhigarhi$dayOfYear, day)
tempData <- weatherData_rakhigarhi[weatherData_rakhigarhi$day == day,]
tempWatbalData <- watbalOutput_rakhigarhi[watbalOutput_rakhigarhi$day == day,]
# solar radiation
weatherSummary_rakhigarhi$solarRadiation.mean <- c(
weatherSummary_rakhigarhi$solarRadiation.mean,
mean(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.sd <- c(
weatherSummary_rakhigarhi$solarRadiation.sd,
sd(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.max <- c(
weatherSummary_rakhigarhi$solarRadiation.max,
max(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.min <- c(
weatherSummary_rakhigarhi$solarRadiation.min,
min(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.error <- c(
weatherSummary_rakhigarhi$solarRadiation.error,
qt(0.975, length(tempData$I) - 1) *
sd(tempData$I, na.rm = TRUE) /
sqrt(length(tempData$I)))
# temperature
## daily mean
weatherSummary_rakhigarhi$temperature.mean <- c(
weatherSummary_rakhigarhi$temperature.mean,
mean(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.sd <- c(
weatherSummary_rakhigarhi$temperature.sd,
sd(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.max <- c(
weatherSummary_rakhigarhi$temperature.max,
max(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.min <- c(
weatherSummary_rakhigarhi$temperature.min,
min(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.error <- c(
weatherSummary_rakhigarhi$temperature.error,
qt(0.975, length(tempData$T2M) - 1) *
sd(tempData$T2M, na.rm = TRUE) /
sqrt(length(tempData$T2M)))
## daily max
weatherSummary_rakhigarhi$maxTemperature.mean <- c(
weatherSummary_rakhigarhi$maxTemperature.mean,
mean(tempData$Tmax, na.rm = TRUE))
weatherSummary_rakhigarhi$maxTemperature.max <- c(
weatherSummary_rakhigarhi$maxTemperature.max,
max(tempData$Tmax, na.rm = TRUE))
weatherSummary_rakhigarhi$maxTemperature.min <- c(
weatherSummary_rakhigarhi$maxTemperature.min,
min(tempData$Tmax, na.rm = TRUE))
weatherSummary_rakhigarhi$maxTemperature.error <- c(
weatherSummary_rakhigarhi$maxTemperature.error,
qt(0.975, length(tempData$Tmax) - 1) *
sd(tempData$Tmax, na.rm = TRUE) /
sqrt(length(tempData$Tmax)))
## daily min
weatherSummary_rakhigarhi$minTemperature.mean <- c(
weatherSummary_rakhigarhi$minTemperature.mean,
mean(tempData$Tmin, na.rm = TRUE))
weatherSummary_rakhigarhi$minTemperature.max <- c(
weatherSummary_rakhigarhi$minTemperature.max,
max(tempData$Tmin, na.rm = TRUE))
weatherSummary_rakhigarhi$minTemperature.min <- c(
weatherSummary_rakhigarhi$minTemperature.min,
min(tempData$Tmin, na.rm = TRUE))
weatherSummary_rakhigarhi$minTemperature.error <- c(
weatherSummary_rakhigarhi$minTemperature.error,
qt(0.975, length(tempData$Tmin) - 1) *
sd(tempData$Tmin, na.rm = TRUE) /
sqrt(length(tempData$Tmin)))
## daily lower and upper deviation
weatherSummary_rakhigarhi$temperature.lowerDeviation <- c(
weatherSummary_rakhigarhi$temperature.lowerDeviation,
mean(tempData$T2M - tempData$Tmin))
weatherSummary_rakhigarhi$temperature.lowerDeviation.error <- c(
weatherSummary_rakhigarhi$temperature.lowerDeviation.error,
qt(0.975, length(tempData$Tmin) - 1) *
sd(tempData$T2M - tempData$Tmin, na.rm = TRUE) /
sqrt(length(tempData$Tmin)))
weatherSummary_rakhigarhi$temperature.upperDeviation <- c(
weatherSummary_rakhigarhi$temperature.upperDeviation,
mean(tempData$Tmax - tempData$T2M))
weatherSummary_rakhigarhi$temperature.upperDeviation.error <- c(
weatherSummary_rakhigarhi$temperature.upperDeviation.error,
qt(0.975, length(tempData$Tmax) - 1) *
sd(tempData$Tmax - tempData$T2M, na.rm = TRUE) /
sqrt(length(tempData$Tmax)))
# precipitation
weatherSummary_rakhigarhi$precipitation.mean <- c(
weatherSummary_rakhigarhi$precipitation.mean,
mean(tempData$RAIN, na.rm = TRUE))
weatherSummary_rakhigarhi$precipitation.max <- c(
weatherSummary_rakhigarhi$precipitation.max,
max(tempData$RAIN, na.rm = TRUE))
weatherSummary_rakhigarhi$precipitation.min <- c(
weatherSummary_rakhigarhi$precipitation.min,
min(tempData$RAIN, na.rm = TRUE))
weatherSummary_rakhigarhi$precipitation.error <- c(
weatherSummary_rakhigarhi$precipitation.error,
qt(0.975, length(tempData$RAIN) - 1) *
sd(tempData$RAIN, na.rm = TRUE) /
sqrt(length(tempData$RAIN)))
# ARID
weatherSummary_rakhigarhi$ARID.mean <- c(
weatherSummary_rakhigarhi$ARID.mean,
mean(tempWatbalData$ARID, na.rm = TRUE))
weatherSummary_rakhigarhi$ARID.max <- c(
weatherSummary_rakhigarhi$ARID.max,
max(tempWatbalData$ARID, na.rm = TRUE))
weatherSummary_rakhigarhi$ARID.min <- c(
weatherSummary_rakhigarhi$ARID.min,
min(tempWatbalData$ARID, na.rm = TRUE))
weatherSummary_rakhigarhi$ARID.error <- c(
weatherSummary_rakhigarhi$ARID.error,
qt(0.975,
length(tempWatbalData$ARID) - 1) *
sd(tempWatbalData$ARID, na.rm = TRUE) /
sqrt(length(tempWatbalData$ARID)))
}
rm(day, tempData, tempWatbalData)
weatherSummary_rakhigarhi <- data.frame(weatherSummary_rakhigarhi)Get ranges of weather variables:
roundToMultiple <- function(i, baseOfMultiple, roundFunction = round)
{
return(match.fun(roundFunction)(i/baseOfMultiple) * baseOfMultiple)
}
rangeSolar = c(
roundToMultiple(min(c(weatherSummary$solarRadiation.min, weatherSummary_rakhigarhi$solarRadiation.min)), 5, floor),
roundToMultiple(max(c(weatherSummary$solarRadiation.max, weatherSummary_rakhigarhi$solarRadiation.max)), 5, ceiling))
rangeTemp = c(
roundToMultiple(min(c(weatherSummary$minTemperature.min, weatherSummary_rakhigarhi$minTemperature.min)), 5, floor),
roundToMultiple(max(c(weatherSummary$maxTemperature.max, weatherSummary_rakhigarhi$maxTemperature.max)), 5, ceiling))
rangePrecip = c(
roundToMultiple(min(c(weatherSummary$precipitation.min, weatherSummary_rakhigarhi$precipitation.min)), 5, floor),
roundToMultiple(max(c(weatherSummary$precipitation.max, weatherSummary_rakhigarhi$precipitation.max)), 5, ceiling))Set colours for maximum and minimum temperature:
maxTemperatureColour = hsv(7.3/360, 74.6/100, 70/100)
minTemperatureColour = hsv(232/360, 64.6/100, 73/100)Define auxiliar function for calculating cumulative curve (extracted from weatherModel.R; see documentation folder in Weather model):
getCumulativePrecipitationOfYear <- function(dailyPrecipitationYear)
{
return(cumsum(dailyPrecipitationYear) / sum(dailyPrecipitationYear))
}Fig - Weather variables
Create figure:
yearLengthInDays_sim = nlevels(factor(weatherSummary$dayOfYear))
yearLengthInDays_real = nlevels(factor(weatherSummary_rakhigarhi$dayOfYear))
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_weatherAndARID.png"
grScale = 2
fontRescale = 0
axisTextRescale = 0
marginTextRescale = 0
leftPlotMargin = 2
rightPlotMargin = 4
cumPrecipMargins <- c(8, 1.5)
bottomMarginWeight = 2
cropRangesHeadLength = 0.1
cropRangesLineWidth = 3
cropRangesLabelsVerticalSpread = 14
cropRangesLabelsSize = 0.95
png(plotName, width = grScale * 800, height = grScale * 1000)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_weatherAndARID.eps"
grScale = 1.2
fontRescale = 0.1
axisTextRescale = -0.1
marginTextRescale = -0.5
leftPlotMargin = 1.5
rightPlotMargin = 3.8
cumPrecipMargins <- c(6, 1.5)
bottomMarginWeight = 1
cropRangesHeadLength = 0.08
cropRangesLineWidth = 3
cropRangesLabelsVerticalSpread = 10
cropRangesLabelsSize = 0.9
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
pointsize = 12,
width = grScale * 8,
height = grScale * 10,
onefile = FALSE,
fallback_resolution = 600,
family = "sans"
)
}
nColumns = 3
nRowsExceptBottom = 5
layout(rbind(matrix(1:(nColumns * nRowsExceptBottom),
nrow = nRowsExceptBottom,
ncol = nColumns,
byrow = FALSE),
c((nColumns * nRowsExceptBottom) + 1,
rep(((nColumns * nRowsExceptBottom) + 2), 2))),
widths = c(1, 10, 12),
heights = c(rep(10, nRowsExceptBottom) + c(1, rep(0, nRowsExceptBottom - 2), -2), bottomMarginWeight)
)
yLabs <- c(expression(paste("solar radiation (", MJ/m^-2, ")")),
"temperature (C)", "precipitation (mm)", "ARID", "crop seasons")
par(cex = grScale,
cex.axis = grScale * (0.8 + axisTextRescale))
# First column: y axis titles
par(mar = c(0, 0, 0, 0.4))
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.5, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[1])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.5, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[2])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.4, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[3])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.6, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[4])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.6, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[5])
#==============================================================================
# Second column: real-world data
# 1. solar radiation
par(mar = c(0.1, leftPlotMargin, 0.1, 0.1))
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$solarRadiation.mean,
axes = FALSE,
ylim = rangeSolar,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$solarRadiation.mean + weatherSummary_rakhigarhi$solarRadiation.error),
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$solarRadiation.mean - weatherSummary_rakhigarhi$solarRadiation.error),
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$solarRadiation.max,
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$solarRadiation.min,
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(rangeSolar[1], rangeSolar[2], 5))
# 2. temperature (daily mean, max, min)
## daily mean
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$temperature.mean,
axes = FALSE,
ylim = rangeTemp,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$temperature.mean + weatherSummary_rakhigarhi$temperature.error),
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$temperature.mean - weatherSummary_rakhigarhi$temperature.error),
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$temperature.max,
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$temperature.min,
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
## daily max
lines(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$maxTemperature.mean,
lwd = grScale, col = maxTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$maxTemperature.mean + weatherSummary_rakhigarhi$maxTemperature.error),
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$maxTemperature.mean - weatherSummary_rakhigarhi$maxTemperature.error),
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$maxTemperature.max,
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$maxTemperature.min,
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
## daily min
lines(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$minTemperature.mean,
lwd = grScale, col = minTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$minTemperature.mean + weatherSummary_rakhigarhi$minTemperature.error),
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$minTemperature.mean - weatherSummary_rakhigarhi$minTemperature.error),
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$minTemperature.max,
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$minTemperature.min,
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(rangeTemp[1], rangeTemp[2], 5))
# 3. precipitation
# cumulative curve
par(mar = c(cumPrecipMargins[1], leftPlotMargin, 0.1, 0.1))
plot(c(1, yearLengthInDays_real), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (year in unique(weatherData_rakhigarhi$year))
{
lines(1:nrow(weatherData_rakhigarhi[weatherData_rakhigarhi$year == year,]),
getCumulativePrecipitationOfYear(weatherData_rakhigarhi[weatherData_rakhigarhi$year == year, "RAIN"]),
lwd = grScale,
col = rgb(0, 0, 0, alpha = 0.05))
}
rm(year)
# daily values
par(new = T,
mar = c(0.1, leftPlotMargin, 0.1, 0.1))
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$precipitation.mean,
axes = FALSE,
ylim = rangePrecip * c(1, cumPrecipMargins[2]),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$precipitation.mean + weatherSummary_rakhigarhi$precipitation.error),
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$precipitation.mean - weatherSummary_rakhigarhi$precipitation.error),
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$precipitation.max,
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$precipitation.min,
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(rangePrecip[1], rangePrecip[2], 10))
# 4. ARID
par(mar = c(0.1, leftPlotMargin, 0.1, 0.1))
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$ARID.mean,
axes = FALSE,
ylim = c(0, 1),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$ARID.mean + weatherSummary_rakhigarhi$ARID.error),
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$ARID.mean - weatherSummary_rakhigarhi$ARID.error),
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$ARID.max,
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$ARID.min,
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(0, 1, 0.1))
#5. crop growing seasons
par(mar = c(3, leftPlotMargin, 0.1, 0.1))
plot(c(1, 365), c(1, nrow(cropTable)), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (aCrop in levels(cropTable$crop))
{
tempData <- cropTable[cropTable$crop == aCrop,]
aCropIndex <- match(aCrop, levels(cropTable$crop))
# winter crop
if (tempData$sugSowingDay > tempData$sugHarvestingDay)
{
arrows(x0 = 1, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = 365, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
else
{
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
}
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# x axis
axis(1, at = c(0,
31,
31+28,
31+28+31,
31+28+31+30,
31+28+31+30+31,
31+28+31+30+31+30,
31+28+31+30+31+30+31,
31+28+31+30+31+30+31+31,
31+28+31+30+31+30+31+31+30,
31+28+31+30+31+30+31+31+30+31,
31+28+31+30+31+30+31+31+30+31+30,
31+28+31+30+31+30+31+31+30+31+30+31),
las = 2
)
#==============================================================================
# Third column: simulated data
# 1. solar radiation
par(mar = c(0.1, 0.1, 0.1, rightPlotMargin))
plot(1:yearLengthInDays_sim,
weatherSummary$solarRadiation.mean,
axes = FALSE,
ylim = rangeSolar,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$solarRadiation.mean + weatherSummary$solarRadiation.error),
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$solarRadiation.mean - weatherSummary$solarRadiation.error),
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$solarRadiation.max,
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$solarRadiation.min,
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# 2. temperature (daily mean, max, min)
## daily mean
plot(1:yearLengthInDays_sim,
weatherSummary$temperature.mean,
axes = FALSE,
ylim = rangeTemp,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$temperature.mean + weatherSummary$temperature.error),
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$temperature.mean - weatherSummary$temperature.error),
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$temperature.max,
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$temperature.min,
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
## daily max
lines(1:yearLengthInDays_sim,
weatherSummary$maxTemperature.mean,
lwd = grScale, col = maxTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$maxTemperature.mean + weatherSummary$maxTemperature.error),
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$maxTemperature.mean - weatherSummary$maxTemperature.error),
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$maxTemperature.max,
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$maxTemperature.min,
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
## daily min
lines(1:yearLengthInDays_sim,
weatherSummary$minTemperature.mean,
lwd = grScale, col = minTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$minTemperature.mean + weatherSummary$minTemperature.error),
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$minTemperature.mean - weatherSummary$minTemperature.error),
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$minTemperature.max,
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$minTemperature.min,
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# 3. precipitation
# cumulative curve
par(mar = c(cumPrecipMargins[1], 0.1, 0.1, rightPlotMargin))
plot(c(1, yearLengthInDays_sim), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (randomSeed in unique(weatherData$randomSeed))
{
for (year in unique(weatherData$currentYear))
{
lines(1:nrow(weatherData[weatherData$randomSeed == randomSeed & weatherData$currentYear == year,]),
getCumulativePrecipitationOfYear(weatherData[weatherData$randomSeed == randomSeed & weatherData$currentYear == year, "RAIN"]),
lwd = grScale,
col = rgb(0, 0, 0, alpha = 0.05))
}
}
rm(year, randomSeed)
axis(4, at = seq(0, 1, 0.25))
mtext("cumulative annual sum", 4, line = 2,
cex = grScale * (1.5 + marginTextRescale))
# daily values
par(new = T,
mar = c(0.1, 0.1, 0.1, rightPlotMargin))
plot(1:yearLengthInDays_sim,
weatherSummary$precipitation.mean,
axes = FALSE,
ylim = rangePrecip * c(1, cumPrecipMargins[2]),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$precipitation.mean + weatherSummary$precipitation.error),
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$precipitation.mean - weatherSummary$precipitation.error),
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$precipitation.max,
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$precipitation.min,
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# 4. ARID
par(mar = c(0.1, 0.1, 0.1, rightPlotMargin))
plot(1:yearLengthInDays_sim,
weatherSummary$ARID.mean,
axes = FALSE,
ylim = c(0, 1),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$ARID.mean + weatherSummary$ARID.error),
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$ARID.mean - weatherSummary$ARID.error),
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$ARID.max,
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$ARID.min,
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
#5. crop growing seasons
par(mar = c(3, 0.1, 0.1, rightPlotMargin))
plot(c(1, 365), c(1, nrow(cropTable)), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (aCrop in levels(cropTable$crop))
{
tempData <- cropTable[cropTable$crop == aCrop,]
aCropIndex <- match(aCrop, levels(cropTable$crop))
# winter crop
if (tempData$sugSowingDay > tempData$sugHarvestingDay)
{
arrows(x0 = 1, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = yearLengthInDays_sim, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
else
{
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
mtext(aCrop,
side = 4, col = cropColours[aCropIndex],
padj = ((1 - aCropIndex/nrow(cropTable)) - 0.4) * cropRangesLabelsVerticalSpread,
las = 2, cex = grScale * cropRangesLabelsSize)
}
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# x axis
axis(1, at = c(0,
31,
31+28,
31+28+31,
31+28+31+30,
31+28+31+30+31,
31+28+31+30+31+30,
31+28+31+30+31+30+31,
31+28+31+30+31+30+31+31,
31+28+31+30+31+30+31+31+30,
31+28+31+30+31+30+31+31+30+31,
31+28+31+30+31+30+31+31+30+31+30,
31+28+31+30+31+30+31+31+30+31+30+31),
las = 2
)
# bottom row: empty and "day of year" or x axis title
par(mar = c(0, 0, 0, 0))
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.45, y = 0.5, font = 4,
cex = grScale * (0.8 + fontRescale),
labels = "day of year")
dev.off()
}
rm(imageFormat, grScale, fontRescale, yLabs, nColumns, nRowsExceptBottom,
axisTextRescale, marginTextRescale, leftPlotMargin, rightPlotMargin, cumPrecipMargins, bottomMarginWeight,
tempData, aCrop, aCropIndex, cropRangesHeadLength, cropRangesLineWidth)knitr::include_graphics(plotName)
Display weather parameters:
knitr::kable(t(yieldData[1, 2:26]), col.names = "default (NW India echo)")| default (NW India echo) | |
|---|---|
| temperature_annualMaxAt2m | 37.00 |
| temperature_annualMinAt2m | 12.80 |
| temperature_meanDailyFluctuation | 2.20 |
| temperature_dailyLowerDeviation | 6.80 |
| temperature_dailyUpperDeviation | 7.90 |
| solar_annualMax | 24.20 |
| solar_annualMin | 9.20 |
| solar_meanDailyFluctuation | 3.30 |
| CO2_annualMin | 245.00 |
| CO2_annualMax | 255.00 |
| CO2_meanDailyFluctuation | 1.00 |
| precipitation_yearlyMean | 489.00 |
| precipitation_yearlySd | 142.20 |
| precipitation_dailyCum_nSamples | 200.00 |
| precipitation_dailyCum_maxSampleSize | 10.00 |
| precipitation_dailyCum_plateauValue_yearlyMean | 0.25 |
| precipitation_dailyCum_plateauValue_yearlySd | 0.10 |
| precipitation_dailyCum_inflection1_yearlyMean | 40.00 |
| precipitation_dailyCum_inflection1_yearlySd | 5.00 |
| precipitation_dailyCum_rate1_yearlyMean | 0.07 |
| precipitation_dailyCum_rate1_yearlySd | 0.02 |
| precipitation_dailyCum_inflection2_yearlyMean | 240.00 |
| precipitation_dailyCum_inflection2_yearlySd | 20.00 |
| precipitation_dailyCum_rate2_yearlyMean | 0.08 |
| precipitation_dailyCum_rate2_yearlySd | 0.02 |
Fig - productivity (yield) per crop type (with default parameter setting)
Summary statistics per crop:
yieldData_summary <- tapply(as.numeric(as.character(yieldData$yield)), yieldData$crop, summary)
yieldData_summary <- data.frame(Reduce(rbind, yieldData_summary), row.names = names(yieldData_summary))## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)
## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)
## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)knitr::kable(yieldData_summary)| Min. | X1st.Qu. | Median | Mean | X3rd.Qu. | Max. | |
|---|---|---|---|---|---|---|
| proso millet | 885.6262 | 1059.3122 | 1109.8132 | 1108.5131 | 1157.0446 | 1373.1612 |
| pearl millet | 300.5779 | 332.6118 | 342.6917 | 342.9984 | 352.8890 | 392.0383 |
| rice | 0.0000 | 230.7686 | 265.1131 | 261.6839 | 299.6243 | 397.6978 |
| barley | 337.8926 | 389.9014 | 412.2040 | 414.7265 | 437.7929 | 520.7963 |
| wheat 1 | 280.4535 | 327.1036 | 345.8900 | 347.8794 | 367.3557 | 437.3639 |
| wheat 2 | 296.1396 | 338.7688 | 357.1001 | 359.2555 | 378.5212 | 447.0028 |
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_yieldPerCrop.png"
grScale = 2
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_yieldPerCrop.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
par(mar = c(6,6,1,1), cex.lab = grScale, cex.axis = 0.8 * grScale)
boxplot(yield ~ factor(crop), data = yieldData,
ylab = expression(paste("yield (", g/m^2, ")")), xlab = "",
col = cropColours
)
dev.off()
}
rm(imageFormat, grScale)knitr::include_graphics(plotName)
ggplot2 option:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_yieldPerCrop_ggplot.png"
grScale = 2
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_yieldPerCrop_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
print({
ggplot(yieldData, aes(y=yield, x=crop, fill=crop)) +
geom_violin(position="dodge") +
#coord_flip() +
#geom_boxplot() +
scale_fill_manual(values = cropColours, name="", guide = FALSE) +
scale_x_discrete() +
theme_light(base_size = 11 * grScale) +
ylab(expression(paste("yield (", g/m^2, ")")))
})
dev.off()
}## Warning: Removed 75 rows containing non-finite values (stat_ydensity).
## Warning: Removed 75 rows containing non-finite values (stat_ydensity).rm(imageFormat, grScale)knitr::include_graphics(plotName)
Fig - mean daily ARID during each crop growing season versus mean ARID of the harvest calendar year
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_ARIDVsARID_grow.png"
grScale = 2
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_ARIDVsARID_grow.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
layout(matrix(c(1, 2), nrow = 1, ncol = 2),
widths = c(10, 3))
par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
boxplot(meanARID_grow ~ factor(crop), data = yieldData,#[yieldData$yield > 0,], # show only non-zero yield
ylab = "mean ARID", xlab = "growing seasons",
#las = 2,
col = cropColours)
text(x = 0.5, y = minARID * 1.005, labels = "a", font = 2, cex = grScale)
par(mar = c(5,1,1,1))
xhist <- hist(yieldData$meanARID, breaks = 20, plot = FALSE)
barplot(xhist$counts, axes = FALSE, space = 0, horiz = TRUE,
xlab = "current year")
text(x = max(xhist$counts) * 0.9, y = minARID * 1.005, labels = "b", font = 2, cex = grScale)
dev.off()
}
rm(imageFormat, grScale)knitr::include_graphics(plotName)
ggplot2 option:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_ARIDVsARID_grow_ggplot.png"
grScale = 2
png(plotName, width = grScale * 600, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_ARIDVsARID_grow_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 6 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
tempData <- rbind(yieldData[, c("crop", "meanARID_grow")],
cbind(crop = rep("year", nrow(yieldData)),
meanARID_grow = yieldData$meanARID))
tempData$meanARID_grow <- as.numeric(tempData$meanARID_grow)
print({
(ggplot(tempData, aes(fill = crop, y = meanARID_grow, x = crop)) +
geom_violin(position = "dodge") +
#geom_boxplot() +
scale_fill_manual(values = c(cropColours, "grey"), name="", guide = FALSE) +
scale_x_discrete() +
theme_light(base_size = 11 * grScale) +
theme(axis.title.x = element_blank()) +
ylab("mean ARID")) +
ylim(minARID, maxARID)
})
dev.off()
}## Warning: Removed 76 rows containing non-finite values (stat_ydensity).
## Warning: Removed 76 rows containing non-finite values (stat_ydensity).rm(tempData, imageFormat, grScale)knitr::include_graphics(plotName)
Fig - mean ARID vs productivity (yield)
Plotting absolute values of yield vs ARID:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_yieldVsARID.png"
grScale = 2
cex.points = 1
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_yieldVsARID.eps"
grScale = 2
cex.points = 0.5
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
layout(matrix(c(1, 2, 3, 3), nrow = 2, ncol = 2, byrow = FALSE),
widths = c(10, 3.5))
par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
plot(c(minARID, maxARID),
c(minYield, maxYield),
xlab = "mean ARID (harvest year)",
ylab = expression(paste("yield (", g/m^2, ")")))
for (aCrop in levels(yieldData$crop))
{
points(yieldData$meanARID[yieldData$crop == aCrop],
yieldData$yield[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(yield ~ meanARID, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minARID * 1.005, y = maxYield * 0.95, labels = "a", font = 2, cex = grScale)
plot(c(minARID, maxARID),
c(minYield, maxYield),
xlab = "mean ARID (grow season)",
ylab = expression(paste("yield (", g/m^2, ")")))
for (aCrop in levels(yieldData$crop))
{
points(yieldData$meanARID_grow[yieldData$crop == aCrop],
yieldData$yield[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(yield ~ meanARID_grow, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minARID * 1.005, y = maxYield * 0.95, labels = "b", font = 2, cex = grScale)
par(mar = c(0, 0, 0, 0), cex = 0.8 * grScale)
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
legend(x = 0, y = 1,
legend = stringi::stri_c(cropTable$S_water[order(match(cropTable$crop, levels(cropTable$crop)))], " (", levels(yieldData$crop), ")"),
title = "S_water",
fill = cropColours)
dev.off()
}
rm(imageFormat, grScale)knitr::include_graphics(plotName)
ggplot alternative:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_yieldVsARID_ggplot.png"
grScale = 2
png(plotName, width = grScale * 350, height = grScale * 500)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_yieldVsARID_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 3.5 * grScale,
height = 5 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
print({
ggplot(data = yieldData, aes(x = meanARID_grow, y = yield, color = crop, fill = crop)) +
geom_point(size = 1) +
geom_smooth(method = 'lm') +
scale_color_manual(values = cropColours, name="", guide = FALSE) +
scale_fill_manual(values = cropColours, name="", guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
xlab("mean ARID") +
ylab(expression(paste("yield (", g/m^2, ")"))) +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))
})
dev.off()
}## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).rm(imageFormat, grScale)knitr::include_graphics(plotName)
Calculate correlation and linear models fitness:
linearModels.table <- data.frame(
levels(yieldData$crop),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop))
)
names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")
for (aCrop in levels(yieldData$crop))
{
linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <-
cor(x = yieldData$meanARID_grow[yieldData$crop == aCrop], y = yieldData$yield[yieldData$crop == aCrop], use = "pairwise.complete.obs")
linearModel <- lm(yield ~ meanARID_grow, data = yieldData[yieldData$crop == aCrop,])
linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop)knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between ARID (x) and Yield (y) per crop")| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | 0.0201669 | 1086.7516 | 0 | 23.79756 | 0.5813393 | -0.0009296 |
| pearl millet | -0.0482022 | 347.6334 | 0 | -7.44716 | 0.1872903 | 0.0009897 |
| rice | -0.9094099 | 693.8981 | 0 | -637.77881 | 0.0000000 | 0.8267950 |
| barley | -0.8629164 | 629.6293 | 0 | -273.61862 | 0.0000000 | 0.7442715 |
| wheat 1 | -0.8470997 | 531.3303 | 0 | -231.94635 | 0.0000000 | 0.7171872 |
| wheat 2 | -0.8582940 | 540.0265 | 0 | -228.55794 | 0.0000000 | 0.7363044 |
Clear temporary objects:
rm(linearModel, linearModels.table,
watbalOutput_rakhigarhi, watbalParameters, soil.FC, soil.WP,
weatherData, weatherData_rakhigarhi, weatherSummary, weatherSummary_rakhigarhi, weatherSummaryNames,
xhist, yieldData, yieldData_summary,
yearLengthInDays_real, yearLengthInDays_sim,
maxARID, minARID, maxYield, minYield, rangePrecip, rangeSolar, rangeTemp,
maxTemperatureColour, minTemperatureColour,
plotName)## Warning in rm(linearModel, linearModels.table, watbalOutput_rakhigarhi, : object
## 'soil.FC' not found## Warning in rm(linearModel, linearModels.table, watbalOutput_rakhigarhi, : object
## 'soil.WP' not foundExperiment 2
Experiment 2 perform the same runs of experiment 1, except by varying stochastically two weather model parameters concerning precipitation: precipitation_yearlyMean (average annual sum of precipitation) and precipitation_dailyCum_plateauValue_yearlyMean (average plateau value of the cumulative curve of daily precipitation).
Loading and preparation
Load data:
yieldData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_type-of-experiment=precipitation-variation_experiment-name=exp2_initRandomSeed=0.csv")Convert yield = 0 to NA when meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):
yieldData$yield[is.na(yieldData$meanARID_grow)] <- NAForce the same order of crops used by cropTable:
yieldData$crop <- factor(yieldData$crop, levels = levels(cropTable$crop))Calculate the range of ARID, yield, precipitation annual sum (yearly mean) and the “plateau value” (yearly mean) of the cumulative curve of daily precipitation:
minARID = round(min(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
maxARID = round(max(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
minYield = round(min(yieldData$yield, na.rm = TRUE), digits = -1)
maxYield = round(max(yieldData$yield, na.rm = TRUE), digits = -1)
minPrecipitation_yearlyMean = round(min(yieldData$precipitation_yearlyMean), digits = -2)
maxPrecipitation_yearlyMean = round(max(yieldData$precipitation_yearlyMean), digits = -2)
minPrecipitation_dailyCum_plateauValue_yearlyMean = round(min(yieldData$precipitation_dailyCum_plateauValue_yearlyMean), digits = 2)
maxPrecipitation_dailyCum_plateauValue_yearlyMean = round(max(yieldData$precipitation_dailyCum_plateauValue_yearlyMean), digits = 2)Display explored values of precipitation_yearlyMean and precipitation_dailyCum_plateauValue_yearlyMean:
layout(matrix(1:2, ncol = 2))
hist(yieldData$precipitation_yearlyMean, main = "precipitation_yearlyMean")
hist(yieldData$precipitation_dailyCum_plateauValue_yearlyMean, main = "precipitation_dailyCum_plateauValue_yearlyMean")
Fig - Effect of precipitation parametres on ARID
Plotting ARID vs precipitation yearly mean and seasonal balance (plateau value):
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp2_ARIDvsPrecipitation.png"
grScale = 2
cex.points = 1
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp2_ARIDvsPrecipitation.eps"
grScale = 2
cex.points = 0.5
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
layout(matrix(c(1, 2, 3, 3), nrow = 2, ncol = 2, byrow = FALSE),
widths = c(10, 3.5))
par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
plot(c(minPrecipitation_yearlyMean, maxPrecipitation_yearlyMean),
c(minARID, maxARID),
xlab = "precipitation annual sum (mm) (yearly mean)",
ylab = "mean ARID (growing seasons)")
for (aCrop in levels(yieldData$crop))
{
points(yieldData$precipitation_yearlyMean[yieldData$crop == aCrop],
yieldData$meanARID_grow[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(meanARID_grow ~ precipitation_yearlyMean, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minPrecipitation_yearlyMean * 0.95, y = minARID * 1.2, labels = "a", font = 2, cex = grScale)
plot(c(minPrecipitation_dailyCum_plateauValue_yearlyMean, maxPrecipitation_dailyCum_plateauValue_yearlyMean),
c(minARID, maxARID),
xlab = "cumulative daily precipitation plateau value (yearly mean)",
ylab = "mean ARID (growing seasons)")
for (aCrop in levels(yieldData$crop))
{
points(yieldData$precipitation_dailyCum_plateauValue_yearlyMean[yieldData$crop == aCrop],
yieldData$meanARID_grow[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(meanARID_grow ~ precipitation_dailyCum_plateauValue_yearlyMean, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minPrecipitation_dailyCum_plateauValue_yearlyMean * 0.95, y = minARID * 1.2, labels = "b", font = 2, cex = grScale)
par(mar = c(0, 0, 0, 0), cex = 0.8 * grScale)
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
legend(x = 0, y = 1,
legend = stringi::stri_c(cropTable$sugSowingDay[order(match(cropTable$crop, levels(cropTable$crop)))], " (", levels(yieldData$crop), ")"),
title = "sowing day of year",
fill = cropColours)
dev.off()
}
rm(imageFormat, grScale, cex.points)knitr::include_graphics(plotName)
ggplot alternative:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp2_ARIDvsPrecipitation_ggplot.png"
grScale = 2
png(plotName, width = grScale * 600, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp2_ARIDvsPrecipitation_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
print({
(ggplot(data = yieldData, aes(y = meanARID_grow, x = precipitation_yearlyMean, color = crop)) +
geom_point(size = 1) +
geom_smooth(method = 'lm') +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2) +
ylab("mean ARID") +
xlab("precipitation annual sum (mm) (yearly mean)") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15),
axis.title = element_text(size = 15),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
(ggplot(data = yieldData, aes(y = meanARID_grow, x = precipitation_dailyCum_plateauValue_yearlyMean, color = crop)) +
geom_point(size = 1) +
geom_smooth(method = 'lm') +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2) +
ylab("") +
xlab("cumulative daily precipitation plateau value (yearly mean)") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15),
axis.title = element_text(size = 15),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
plot_layout(nrow = 1, widths = c(1, 1))
})
dev.off()
}## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).rm(imageFormat, grScale)knitr::include_graphics(plotName)
Calculate correlation and linear models fitness:
linearModels.table <- data.frame(
levels(yieldData$crop),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop))
)
names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")
for (aCrop in levels(yieldData$crop))
{
linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <-
cor(x = yieldData$meanARID_grow[yieldData$crop == aCrop],
y = yieldData$precipitation_yearlyMean[yieldData$crop == aCrop],
use = "pairwise.complete.obs")
linearModel <- lm(meanARID_grow ~ precipitation_yearlyMean, data = yieldData[yieldData$crop == aCrop,])
linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop)knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between precipitation annual sum yearly mean (x) and ARID (y) per crop")| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | -0.5760686 | 0.9960945 | 0 | -0.0001314 | 0 | 0.3309618 |
| pearl millet | -0.7070320 | 0.9428775 | 0 | -0.0003964 | 0 | 0.4992257 |
| rice | -0.7308760 | 0.9545028 | 0 | -0.0003454 | 0 | 0.5335569 |
| barley | -0.7133386 | 0.8883206 | 0 | -0.0004509 | 0 | 0.5081726 |
| wheat 1 | -0.7138806 | 0.8916387 | 0 | -0.0004425 | 0 | 0.5089473 |
| wheat 2 | -0.7138806 | 0.8916387 | 0 | -0.0004425 | 0 | 0.5089473 |
linearModels.table <- data.frame(
levels(yieldData$crop),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop))
)
names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")
for (aCrop in levels(yieldData$crop))
{
linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <-
cor(x = yieldData$meanARID_grow[yieldData$crop == aCrop],
y = yieldData$precipitation_dailyCum_plateauValue_yearlyMean[yieldData$crop == aCrop],
use = "pairwise.complete.obs")
linearModel <- lm(meanARID_grow ~ precipitation_dailyCum_plateauValue_yearlyMean, data = yieldData[yieldData$crop == aCrop,])
linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop)knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between cumulative daily precipitation plateau value yearly mean (x) and ARID (y) per crop")| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | 0.1633846 | 0.8971123 | 0 | 0.0628231 | 6.9e-06 | 0.0253933 |
| pearl millet | 0.4151101 | 0.5367639 | 0 | 0.3922293 | 0.0e+00 | 0.1712098 |
| rice | 0.4268577 | 0.6015239 | 0 | 0.3400367 | 0.0e+00 | 0.1811142 |
| barley | -0.3736394 | 0.8743004 | 0 | -0.3980747 | 0.0e+00 | 0.1384163 |
| wheat 1 | -0.3754886 | 0.8787469 | 0 | -0.3923104 | 0.0e+00 | 0.1398036 |
| wheat 2 | -0.3754886 | 0.8787469 | 0 | -0.3923104 | 0.0e+00 | 0.1398036 |
Clear temporary objects:
rm(linearModel, linearModels.table,
yieldData,
maxARID, minARID, maxYield, minYield,
maxPrecipitation_yearlyMean, minPrecipitation_yearlyMean,
maxPrecipitation_dailyCum_plateauValue_yearlyMean, minPrecipitation_dailyCum_plateauValue_yearlyMean,
plotName)Fig - terrains generated with the Land model
Loading and preparation
Defining terrain random number generator seed (must be available inside “terrains” folder:
listOfTerrainSeeds <- c("0", "35", "56", "72", "92")Setting dimensions:
figScale = 2
terrainDim = 450 * figScale
topHeight = 50 * figScale
bottomHeight = 150 * figScale
leftWidth = 100 * figScale
internalMargin = 5 * figScale
figWidth = terrainDim * 3 + leftWidth + 2 * internalMargin
figHeight = terrainDim * length(listOfTerrainSeeds) + topHeight + bottomHeight + (length(listOfTerrainSeeds) - 1) * internalMarginCreate separate component images
Create margin header images:
headers <- c("seed",
"Elevation",
"Soil texture",
"Soil texture type")
namesInFile <- c("seed", "elev", "soilText", "soilTextType")
for (i in 1:length(headers))
{
#if (i == 1) { width = leftWidth } else { width = terrainDim }
width = ifelse(i == 1, leftWidth, terrainDim)
png(paste0("FigTerrains_tempFiles/FigTerrains_", namesInFile[i], ".PNG"),
width = width,
height = topHeight)
par(mar = rep(0, 4))
plot(c(0, 1), c(0, 1), type = "n", axes = FALSE, ann = FALSE)
text(0.5, 0.7, cex = figScale * 3.7, adj = c(0.5, 0.5),
labels = headers[i])
dev.off()
}
rm(i)
for (seed in listOfTerrainSeeds)
{
png(paste0("FigTerrains_tempFiles/FigTerrains_seed=", seed, ".PNG"),
width = leftWidth,
height = terrainDim)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
text(0.5, 0.5, cex = figScale * 5, adj = .5,
labels = seed)
dev.off()
}
rm(seed)Create legend images:
# code converted from NetLogo (maximum = 255)
elev <- 100/255 + (155/255) * seq(0, 1, length.out = 10)
png("FigTerrains_tempFiles/FigTerrains_legend_elev.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = c(2,2,1,2), cex = figScale * 2)
barplot(rep(100, 10),
col = rgb(elev - (100/255), elev, 0),
border = rgb(elev - (100/255), elev, 0),
axes = FALSE, space = 0)
mtext("min.", side = 1, line = 0.1, adj = 0, cex = figScale * 2)
mtext("max.", side = 1, line = 0.1, adj = 1, cex = figScale * 2)
#mtext("elevation", side = 1, line = 1.2, adj = .5, cex = figScale * 3)
dev.off()## png
## 2#======================================================================
soilTexture <- c("100% sand", "100% silt", "100% clay")
soilTexture_colors <- c("red", "green", "blue")
xy <- list(c(0.2, 1.05),
c(0.2, 0.75),
c(0.2, 0.45))
png("FigTerrains_tempFiles/FigTerrains_legend_soilText.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
for (i in 1:length(soilTexture))
{
legend(xy[[i]][1], xy[[i]][2],
legend = soilTexture[i],
fill = soilTexture_colors[i],
bty = "n", adj = 0.1,
cex = figScale * 2.5, pt.cex = figScale * 5)
}
rm(i)
dev.off()## png
## 2#======================================================================
soilTextureTypes <- c("Sand", "Loamy sand", "Sandy loam",
"Loam", "Silt loam", "Silty clay",
"Clay", "Clay loam", "Sandy clay loam")
soilTextureTypes_colors <- c("#d73229", "#f16a15", "#9d6e48",
"#eded31", "#59b03c", "#54c4c4",
"#2d8dbe", "#345da9", "#a71b6a")
xy <- list(c(-0.08, 1), c(0.17, 1), c(0.55, 1),
c(-0.08, 0.7), c(0.17, 0.7), c(0.55, 0.7),
c(-0.08, 0.4), c(0.17, 0.4), c(0.55, 0.4))
png("FigTerrains_tempFiles/FigTerrains_legend_soilTextType.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
for (i in 1:length(soilTextureTypes))
{
legend(xy[[i]][1], xy[[i]][2],
legend = soilTextureTypes[i],
fill = soilTextureTypes_colors[i],
bty = "n", adj = 0.1,
cex = figScale * 1.8, pt.cex = figScale * 5)
}
rm(i)
dev.off()## png
## 2Create blank image:
png("FigTerrains_tempFiles/FigTerrains_blank.PNG",
width = 10,
height = 10)## Warning in png("FigTerrains_tempFiles/FigTerrains_blank.PNG", width = 10, :
## 'width=10, height=10' are unlikely values in pixelsdev.off()## png
## 2Compose final figure
Load input images:
elevTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfTerrainSeeds ,".PNG")
soilTextTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfTerrainSeeds ,"_soilTexture.PNG")
soilTextTypeTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfTerrainSeeds ,"_soilTextureTypes.PNG")
numberOfColumns = 6
numberOfRows = 11
# it is important that these plots were saved as "PNG not "png"
listOfPlots <- c("FigTerrains_tempFiles/FigTerrains_seed.PNG",
"FigTerrains_tempFiles/FigTerrains_elev.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_soilText.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_soilTextType.PNG",
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[1], ".PNG"),
elevTerrains[1],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[1],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[1],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[2], ".PNG"),
elevTerrains[2],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[2],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[2],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[3], ".PNG"),
elevTerrains[3],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[3],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[3],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[4], ".PNG"),
elevTerrains[4],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[4],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[4],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[5], ".PNG"),
elevTerrains[5],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[5],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[5],
# ---
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_legend_elev.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_legend_soilText.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_legend_soilTextType.PNG"
)
grobs <- lapply(lapply(listOfPlots, png::readPNG), grid::rasterGrob)Create Figure:
widths = 50 * c(leftWidth/figWidth,
terrainDim/figWidth,
internalMargin/figWidth,
terrainDim/figWidth,
internalMargin/figWidth,
terrainDim/figWidth) * figScale
heights = (figHeight/figWidth) * 50 * c(topHeight/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
bottomHeight/figHeight) * figScale
lay <- matrix(1:(numberOfColumns * numberOfRows), ncol = numberOfColumns, byrow = TRUE)
png("figures/Fig_terrains.png", width = figWidth, height = figHeight)
gridExtra::grid.arrange(grobs = grobs,
layout_matrix = lay,
widths = unit(widths, "cm"),
heights = unit(heights, "cm"))
dev.off()## png
## 2knitr::include_graphics("figures/Fig_terrains.png")